Chapter 7 Functional differences

load("data/data.Rdata")

7.1 Data preparation

# Aggregate bundle-level GIFTs into the compound level
GIFTs_elements <- to.elements(genome_gifts, GIFT_db3)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts$genome, ]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
  select_if(~ !is.numeric(.) || sum(.) != 0)

elements <- GIFTs_elements_filtered %>%
  as.data.frame()

# Aggregate element-level GIFTs into the function level
GIFTs_functions <- to.functions(GIFTs_elements_filtered, GIFT_db3)
functions <- GIFTs_functions %>%
  as.data.frame()

# Aggregate function-level GIFTs into overall Biosynthesis, Degradation and Structural GIFTs
GIFTs_domains <- to.domains(GIFTs_functions, GIFT_db3)
domains <- GIFTs_domains %>%
  as.data.frame()

# Get community-weighed average GIFTs per sample
GIFTs_elements_community <- to.community(GIFTs_elements_filtered, genome_counts_filt %>% column_to_rownames(., "genome") %>% tss(), GIFT_db3)
GIFTs_functions_community <- to.community(GIFTs_functions, genome_counts_filt %>% column_to_rownames(., "genome") %>% tss(), GIFT_db3)
GIFTs_domains_community <- to.community(GIFTs_domains, genome_counts_filt %>% column_to_rownames(., "genome") %>% tss(), GIFT_db3)

uniqueGIFT_db<- unique(GIFT_db3[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

7.2 Genomes GIFT profiles

GIFTs_elements %>%
  as_tibble(., rownames = "MAG") %>%
  reshape2::melt() %>%
  rename(Code_element = variable, GIFT = value) %>%
  inner_join(GIFT_db3,by="Code_element") %>%
  ggplot(., aes(x=Code_element, y=MAG, fill=GIFT, group=Code_function))+
    geom_tile()+
    scale_y_discrete(guide = guide_axis(check.overlap = TRUE))+
    scale_x_discrete(guide = guide_axis(check.overlap = TRUE))+
    scale_fill_gradientn(colours=rev(c("#d53e4f", "#f46d43", "#fdae61", "#fee08b", "#e6f598", "#abdda4", "#ddf1da")))+
    facet_grid(. ~ Code_function, scales = "free", space = "free")+
    theme_grey(base_size=8)+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),strip.text.x = element_text(angle = 90))

7.3 Function level

GIFTs_functions_community %>%
    as.data.frame() %>%
    rownames_to_column(var="sample") %>%
    filter(sample!="AD69") %>%
    pivot_longer(!sample,names_to="trait",values_to="gift") %>%
    left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
    ggplot(aes(x=trait,y=time_point,fill=gift)) +
        geom_tile(colour="white", size=0.2)+
        scale_fill_gradientn(colours=rev(c("#d53e4f", "#f46d43", "#fdae61", "#fee08b", "#e6f598", "#abdda4", "#ddf1da")))+
        facet_grid(type ~ ., scales="free",space="free")

7.4 Element level

GIFTs_elements_community_merged<-GIFTs_elements_community %>%
    as.data.frame() %>%
    rownames_to_column(var="sample") %>%
    filter(sample!="AD69") %>%
    pivot_longer(!sample,names_to="trait",values_to="gift") %>%
    left_join(sample_metadata, by = join_by(sample == Tube_code))%>%
    mutate(functionid = substr(trait, 1, 3)) %>%
    mutate(trait = case_when(
      trait %in% GIFT_db3$Code_element ~ GIFT_db3$Element[match(trait, GIFT_db3$Code_element)],
      TRUE ~ trait
    )) %>%
    mutate(functionid = case_when(
      functionid %in% GIFT_db3$Code_function ~ GIFT_db3$Function[match(functionid, GIFT_db3$Code_function)],
      TRUE ~ functionid
    )) %>%
    mutate(trait=factor(trait,levels=unique(GIFT_db3$Element))) %>%
    mutate(functionid=factor(functionid,levels=unique(GIFT_db3$Function)))

# Create an interaction variable for time_point and sample
GIFTs_elements_community_merged$interaction_var <- interaction(GIFTs_elements_community_merged$sample, GIFTs_elements_community_merged$time_point)
  
ggplot(GIFTs_elements_community_merged,aes(x=interaction_var,y=trait,fill=gift)) +
        geom_tile(colour="white", linewidth=0.2)+
        scale_fill_gradientn(colours=rev(c("#d53e4f", "#f46d43", "#fdae61", "#fee08b", "#e6f598", "#abdda4", "#ddf1da")))+
        facet_grid(functionid ~ type, scales="free",space="free") +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=5),
              strip.text.y = element_text(angle = 0)) + 
        labs(y="Traits",x="Time_point",fill="GIFT")+
  scale_x_discrete(labels = function(x) gsub(".*\\.", "", x))

7.5 Comparison of samples from the 0 Time_point (0_Wild)

7.5.1 GIFTs Functional community

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  group_by(species) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  species             MCI     sd
  <chr>             <dbl>  <dbl>
1 Podarcis_liolepis 0.271 0.0203
2 Podarcis_muralis  0.281 0.0146

7.5.1.1 GIFT test visualisation

GIFTs_functions_community %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  select(c(1:30, 33)) %>%
  pivot_longer(-c(sample,Population),names_to = "trait", values_to = "value") %>%
  mutate(trait = case_when(
      trait %in% GIFT_db3$Code_function ~ GIFT_db3$Function[match(trait, GIFT_db3$Code_function)],
      TRUE ~ trait
    )) %>%
  mutate(trait=factor(trait,levels=unique(GIFT_db3$Function))) %>%
  ggplot(aes(x=value, y=Population, group=Population, fill=Population, color=Population)) +
    geom_boxplot() +
    scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
    facet_grid(trait ~ ., space="free", scales="free") +
              theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
              strip.text.y = element_text(angle = 0)) + 
        labs(y="Traits",x="Metabolic capacity index")

7.5.2 GIFTs Domain community

GIFTs_domains_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  group_by(species) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  species             MCI     sd
  <chr>             <dbl>  <dbl>
1 Podarcis_liolepis 0.283 0.0162
2 Podarcis_muralis  0.289 0.0184

7.5.3 GIFTs Elements community

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  group_by(species) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  species             MCI     sd
  <chr>             <dbl>  <dbl>
1 Podarcis_liolepis 0.241 0.0239
2 Podarcis_muralis  0.260 0.0163
sample_metadata_wild <- sample_metadata%>% 
  filter(time_point == "0_Wild")

element_gift_wild <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(., sample_metadata_wild[c(1,3)], by="Tube_code")
# Find numeric columns
numeric_cols <- sapply(element_gift_wild, is.numeric)

# Calculate column sums for numeric columns only
col_sums_numeric <- colSums(element_gift_wild[, numeric_cols])

# Identify numeric columns with sums not equal to zero
nonzero_numeric_cols <- names(col_sums_numeric)[col_sums_numeric != 0]

# Remove numeric columns with sums not equal to zero
filtered_data <- element_gift_wild[, !numeric_cols | colnames(element_gift_wild) %in% nonzero_numeric_cols]
significant_elements_wild <- filtered_data %>%
  pivot_longer(-c(Tube_code,species), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ species, exact=FALSE)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift_wild  %>% 
  dplyr::select(-c(species))  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements_wild$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_wild[c(1,3)], by = join_by(Tube_code == Tube_code))

element_gift_filt %>%
  dplyr::select(-Tube_code)%>%
  group_by(species)  %>%
  summarise(across(everything(), mean))%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))

element_gift_names <- element_gift_filt%>%
  dplyr::select(-species)%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
  dplyr::select(-Elements)%>%
  dplyr::select(Function, everything())%>%
  t()%>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_wild[c(1,4)], by = join_by(Tube_code == Tube_code))
colNames <- names(element_gift_names)[2:72] #always check names(element_gift_names) first to know where your traits finish
for(i in colNames){
  plt <- ggplot(element_gift_names, aes(x=Population, y=.data[[i]], color = Population)) +
    geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
    geom_jitter(width = 0.1, show.legend = TRUE) +
    theme_minimal() +
    theme(
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.border = element_blank())
  print(plt)
}

7.6 Comparison of samples from the 1st Time_point (1_Acclimation)

7.6.1 GIFTs Functional community

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  group_by(species) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  species             MCI     sd
  <chr>             <dbl>  <dbl>
1 Podarcis_liolepis 0.284 0.0148
2 Podarcis_muralis  0.273 0.0258

7.6.1.1 GIFT test visualisation

GIFTs_functions_community %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  select(c(1:30, 33)) %>%
  pivot_longer(-c(sample,Population),names_to = "trait", values_to = "value") %>%
  mutate(trait = case_when(
      trait %in% GIFT_db3$Code_function ~ GIFT_db3$Function[match(trait, GIFT_db3$Code_function)],
      TRUE ~ trait
    )) %>%
  mutate(trait=factor(trait,levels=unique(GIFT_db3$Function))) %>%
  ggplot(aes(x=value, y=Population, group=Population, fill=Population, color=Population)) +
    geom_boxplot() +
    scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
    facet_grid(trait ~ ., space="free", scales="free") +
              theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
              strip.text.y = element_text(angle = 0)) + 
        labs(y="Traits",x="Metabolic capacity index")

7.6.2 GIFTs Domain community

GIFTs_domains_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  group_by(species) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  species             MCI     sd
  <chr>             <dbl>  <dbl>
1 Podarcis_liolepis 0.286 0.0130
2 Podarcis_muralis  0.277 0.0306

7.6.3 GIFTs Elements community

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  group_by(species) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 2 × 3
  species             MCI     sd
  <chr>             <dbl>  <dbl>
1 Podarcis_liolepis 0.264 0.0182
2 Podarcis_muralis  0.254 0.0271
sample_metadata_accli <- sample_metadata%>% 
  filter(time_point == "1_Acclimation")

element_gift_accli <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(., sample_metadata_accli[c(1,3)], by="Tube_code")
# Find numeric columns
numeric_cols <- sapply(element_gift_accli, is.numeric)

# Calculate column sums for numeric columns only
col_sums_numeric <- colSums(element_gift_accli[, numeric_cols])

# Identify numeric columns with sums not equal to zero
nonzero_numeric_cols <- names(col_sums_numeric)[col_sums_numeric != 0]

# Remove numeric columns with sums not equal to zero
filtered_data <- element_gift_accli[, !numeric_cols | colnames(element_gift_accli) %in% nonzero_numeric_cols]
significant_elements_accli <- filtered_data %>%
  pivot_longer(-c(Tube_code,species), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ species, exact=FALSE)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift_accli  %>% 
  dplyr::select(-c(species))  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements_accli$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_accli[c(1,3)], by = join_by(Tube_code == Tube_code))

element_gift_filt %>%
  dplyr::select(-Tube_code)%>%
  group_by(species)  %>%
  summarise(across(everything(), mean))%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))

element_gift_names <- element_gift_filt%>%
  dplyr::select(-species)%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
  dplyr::select(-Elements)%>%
  dplyr::select(Function, everything())%>%
  t()%>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_accli[c(1,4)], by = join_by(Tube_code == Tube_code))
colNames <- names(element_gift_names)[2:51] #always check names(element_gift_names) first to know where your traits finish
for(i in colNames){
  plt <- ggplot(element_gift_names, aes(x=Population, y=.data[[i]], color = Population)) +
    geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
    geom_jitter(width = 0.1, show.legend = TRUE) +
    theme_minimal() +
    theme(
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.border = element_blank())
  print(plt)
}

7.7 Comparison of samples from the 5th Time_point (5_Post-FMT1)

7.7.1 GIFTs Functional community

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  group_by(type) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 3 × 3
  type          MCI     sd
  <chr>       <dbl>  <dbl>
1 Control     0.308 0.0247
2 Hot_control 0.309 0.0464
3 Treatment   0.291 0.0120

7.7.1.1 GIFT test visualisation

GIFTs_functions_community %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  select(c(1:30, 36)) %>%
  pivot_longer(-c(sample,type),names_to = "trait", values_to = "value") %>%
  mutate(trait = case_when(
      trait %in% GIFT_db3$Code_function ~ GIFT_db3$Function[match(trait, GIFT_db3$Code_function)],
      TRUE ~ trait
    )) %>%
  mutate(trait=factor(trait,levels=unique(GIFT_db3$Function))) %>%
  ggplot(aes(x=value, y=type, group=type, fill=type, color=type)) +
    geom_boxplot() +
    scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c","#76b183")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
    facet_grid(trait ~ ., space="free", scales="free") +
              theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
              strip.text.y = element_text(angle = 0)) + 
        labs(y="Traits",x="Metabolic capacity index")

7.7.2 GIFTs Domain community

GIFTs_domains_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  group_by(type) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 3 × 3
  type          MCI     sd
  <chr>       <dbl>  <dbl>
1 Control     0.305 0.0230
2 Hot_control 0.307 0.0447
3 Treatment   0.287 0.0141

7.7.3 GIFTs Elements community

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  group_by(type) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 3 × 3
  type          MCI     sd
  <chr>       <dbl>  <dbl>
1 Control     0.293 0.0303
2 Hot_control 0.293 0.0467
3 Treatment   0.275 0.0137
sample_metadata_tm5 <- sample_metadata%>% 
  filter(time_point == "5_Post-FMT1")%>% 
  filter(type != "Hot_control")

element_gift_tm5 <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(sample_metadata_tm5 %>% select(1, 7), by = "Tube_code")
# Find numeric columns
numeric_cols <- sapply(element_gift_tm5, is.numeric)

# Calculate column sums for numeric columns only
col_sums_numeric <- colSums(element_gift_tm5[, numeric_cols])

# Identify numeric columns with sums not equal to zero
nonzero_numeric_cols <- names(col_sums_numeric)[col_sums_numeric != 0]

# Remove numeric columns with sums not equal to zero
filtered_data <- element_gift_tm5[, !numeric_cols | colnames(element_gift_tm5) %in% nonzero_numeric_cols]
significant_elements_tm5 <- filtered_data %>%
  pivot_longer(-c(Tube_code,type), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ type, exact=FALSE)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_value < 0.05)  %>% #take into account that p_value is used and not p_adjust
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift_tm5  %>% 
  dplyr::select(-c(type))  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements_tm5$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_tm5[c(1,7)], by = join_by(Tube_code == Tube_code))

element_gift_filt %>%
  dplyr::select(-Tube_code)%>%
  group_by(type)  %>%
  summarise(across(everything(), mean))%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))

element_gift_names <- element_gift_filt%>%
  dplyr::select(-type)%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
  dplyr::select(-Elements)%>%
  dplyr::select(Function, everything())%>%
  t()%>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_tm5[c(1,7)], by = join_by(Tube_code == Tube_code))
colNames <- names(element_gift_names)[2:14] #always check names(element_gift_names) first to now where your traits finish
for(i in colNames){
  plt <- ggplot(element_gift_names, aes(x=type, y=.data[[i]], color = type)) +
    geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
    geom_jitter(width = 0.1, show.legend = TRUE) +
    theme_minimal() +
    theme(
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.border = element_blank())
  print(plt)
}

7.8 Comparison of samples from the 6th Time_point (6_Post-FMT2)

7.8.1 GIFTs Functional community

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  group_by(type) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 3 × 3
  type          MCI     sd
  <chr>       <dbl>  <dbl>
1 Control     0.287 0.0157
2 Hot_control 0.285 0.0294
3 Treatment   0.280 0.0197
MCI_tm6 <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2")

shapiro.test(MCI_tm6$value)#normality test

    Shapiro-Wilk normality test

data:  MCI_tm6$value
W = 0.95761, p-value = 0.3256
res.aov<-aov(value ~ type, data=MCI_tm6)#anova test
summary(res.aov)
            Df   Sum Sq   Mean Sq F value Pr(>F)
type         2 0.000195 0.0000974   0.195  0.824
Residuals   24 0.011994 0.0004998               

7.8.1.1 GIFT test visualisation

GIFTs_functions_community %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  select(c(1:30, 36)) %>%
  pivot_longer(-c(sample,type),names_to = "trait", values_to = "value") %>%
  mutate(trait = case_when(
      trait %in% GIFT_db3$Code_function ~ GIFT_db3$Function[match(trait, GIFT_db3$Code_function)],
      TRUE ~ trait
    )) %>%
  mutate(trait=factor(trait,levels=unique(GIFT_db3$Function))) %>%
  ggplot(aes(x=value, y=type, group=type, fill=type, color=type)) +
    geom_boxplot() +
    scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c","#76b183")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
    facet_grid(trait ~ ., space="free", scales="free") +
              theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
              strip.text.y = element_text(angle = 0)) + 
        labs(y="Traits",x="Metabolic capacity index")

7.8.2 GIFTs Domain community

GIFTs_domains_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  group_by(type) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 3 × 3
  type          MCI     sd
  <chr>       <dbl>  <dbl>
1 Control     0.290 0.0210
2 Hot_control 0.292 0.0310
3 Treatment   0.282 0.0190
MCI_domains_tm6 <- GIFTs_domains_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2")

shapiro.test(MCI_domains_tm6$value)#normality test

    Shapiro-Wilk normality test

data:  MCI_domains_tm6$value
W = 0.9498, p-value = 0.2119
res.aov<-aov(value ~ type, data=MCI_domains_tm6)#anova test
summary(res.aov)
            Df  Sum Sq   Mean Sq F value Pr(>F)
type         2 0.00057 0.0002851   0.485  0.622
Residuals   24 0.01411 0.0005878               

7.8.3 GIFTs Elements community

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  group_by(type) %>%
  summarise(MCI = mean(value), sd = sd(value))
# A tibble: 3 × 3
  type          MCI     sd
  <chr>       <dbl>  <dbl>
1 Control     0.268 0.0169
2 Hot_control 0.263 0.0314
3 Treatment   0.261 0.0219
MCI_ele_tm6 <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2")

shapiro.test(MCI_ele_tm6$value) #normality test

    Shapiro-Wilk normality test

data:  MCI_ele_tm6$value
W = 0.9806, p-value = 0.8757
res.aov<-aov(value ~ type, data=MCI_ele_tm6) #anova test
summary(res.aov)
            Df   Sum Sq   Mean Sq F value Pr(>F)
type         2 0.000291 0.0001453   0.249  0.782
Residuals   24 0.014001 0.0005834               

7.8.3.1 CI vs CC

sample_metadata_TM6 <- sample_metadata%>% 
  filter(time_point == "6_Post-FMT2")%>% 
  filter(type != "Hot_control")

element_gift_TM6 <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(sample_metadata_TM6 %>% select(1, 7), by = "Tube_code")
# Find numeric columns
numeric_cols <- sapply(element_gift_TM6, is.numeric)

# Calculate column sums for numeric columns only
col_sums_numeric <- colSums(element_gift_TM6[, numeric_cols])

# Identify numeric columns with sums not equal to zero
nonzero_numeric_cols <- names(col_sums_numeric)[col_sums_numeric != 0]

# Remove numeric columns with sums not equal to zero
filtered_data <- element_gift_TM6[, !numeric_cols | colnames(element_gift_TM6) %in% nonzero_numeric_cols]
significant_elements_TM6 <- filtered_data %>%
  pivot_longer(-c(Tube_code,type), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ type, exact=FALSE)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_value < 0.05)  %>% #take into account that p_value is used and not p_adjust
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift_TM6  %>% 
  dplyr::select(-c(type))  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements_TM6$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_TM6[c(1,7)], by = join_by(Tube_code == Tube_code))

element_gift_filt %>%
  dplyr::select(-Tube_code)%>%
  group_by(type)  %>%
  summarise(across(everything(), mean))%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))

element_gift_names <- element_gift_filt%>%
  dplyr::select(-type)%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
  dplyr::select(-Elements)%>%
  dplyr::select(Function, everything())%>%
  t()%>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_TM6[c(1,7)], by = join_by(Tube_code == Tube_code))
colNames <- names(element_gift_names)[2:20] #always check names(element_gift_names) first to now where your traits finish
for(i in colNames){
  plt <- ggplot(element_gift_names, aes(x=type, y=.data[[i]], color = type)) +
    geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
    geom_jitter(width = 0.1, show.legend = TRUE) +
    theme_minimal() +
    theme(
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.border = element_blank())
  print(plt)
}

7.8.3.2 CI vs WC

sample_metadata_TM8 <- sample_metadata%>% 
  filter(time_point == "6_Post-FMT2")%>% 
  filter(type != "Control")

element_gift_TM8 <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(sample_metadata_TM8 %>% select(1, 7), by = "Tube_code")
# Find numeric columns
numeric_cols <- sapply(element_gift_TM8, is.numeric)

# Calculate column sums for numeric columns only
col_sums_numeric <- colSums(element_gift_TM8[, numeric_cols])

# Identify numeric columns with sums not equal to zero
nonzero_numeric_cols <- names(col_sums_numeric)[col_sums_numeric != 0]

# Remove numeric columns with sums not equal to zero
filtered_data <- element_gift_TM8[, !numeric_cols | colnames(element_gift_TM8) %in% nonzero_numeric_cols]
significant_elements_TM8 <- filtered_data %>%
  pivot_longer(-c(Tube_code,type), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ type, exact=FALSE)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_value < 0.05)  %>% #take into account that p_value is used and not p_adjust
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift_TM8  %>% 
  dplyr::select(-c(type))  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements_TM8$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_TM8[c(1,7)], by = join_by(Tube_code == Tube_code))

element_gift_filt %>%
  dplyr::select(-Tube_code)%>%
  group_by(type)  %>%
  summarise(across(everything(), mean))%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))

element_gift_names <- element_gift_filt%>%
  dplyr::select(-type)%>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Elements")  %>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element))%>%
  dplyr::select(-Elements)%>%
  dplyr::select(Function, everything())%>%
  t()%>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "Tube_code")%>% 
  left_join(., sample_metadata_TM8[c(1,7)], by = join_by(Tube_code == Tube_code))
colNames <- names(element_gift_names)[2:13] #always check names(element_gift_names) first to now where your traits finish
for(i in colNames){
  plt <- ggplot(element_gift_names, aes(x=type, y=.data[[i]], color = type)) +
    geom_boxplot(alpha = 0.2, outlier.shape = NA, width = 0.3, show.legend = FALSE) +
    geom_jitter(width = 0.1, show.legend = TRUE) +
    theme_minimal() +
    theme(
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.border = element_blank())
  print(plt)
}

7.9 Domain level

7.9.1 Comparison of samples from the 0 Time_point (0_Wild)

#Merge the functional domains with the metadata
merge_gift_wild<- GIFTs_domains_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(., sample_metadata_wild, by="Tube_code")
#Biosynthesis
p1 <-merge_gift_wild %>%
  ggplot(aes(x=species,y=Biosynthesis,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Degradation
p2 <-merge_gift_wild %>%
  ggplot(aes(x=species,y=Degradation,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.45, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Structure
p3 <-merge_gift_wild %>%
  ggplot(aes(x=species,y=Structure,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Transport
p4 <-merge_gift_wild %>%
  ggplot(aes(x=species,y=Transport,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

7.9.2 Comparison of samples from the 1st Time_point (1_Acclimation)

#Merge the functional domains with the metadata
merge_gift_accli<- GIFTs_domains_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(., sample_metadata_accli, by="Tube_code")
#Biosynthesis
p1 <-merge_gift_accli %>%
  ggplot(aes(x=species,y=Biosynthesis,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Degradation
p2 <-merge_gift_accli %>%
  ggplot(aes(x=species,y=Degradation,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.45, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Structure
p3 <-merge_gift_accli %>%
  ggplot(aes(x=species,y=Structure,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Transport
p4 <-merge_gift_accli %>%
  ggplot(aes(x=species,y=Transport,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

7.9.3 Comparison of samples from the 5th Time_point (5_Post-FMT1)

#Merge the functional domains with the metadata
merge_gift_tm5<- GIFTs_domains_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(., sample_metadata_tm5, by="Tube_code")
#Biosynthesis
p1 <-merge_gift_tm5 %>%
  ggplot(aes(x=type,y=Biosynthesis,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Degradation
p2 <-merge_gift_tm5 %>%
  ggplot(aes(x=type,y=Degradation,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.45, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Structure
p3 <-merge_gift_tm5 %>%
  ggplot(aes(x=type,y=Structure,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

#Transport
p4 <-merge_gift_tm5 %>%
  ggplot(aes(x=species,y=Transport,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")
Warning: Computation failed in `stat_compare_means()`.
Caused by error:
! argument "x" is missing, with no default

7.9.4 Comparison of samples from the 6th Time_point (6_Post-FMT2)

7.9.4.1 CI vs CC

#Merge the functional domains with the metadata
merge_gift_TM6 <- GIFTs_domains_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(., sample_metadata_TM6, by="Tube_code")
#Biosynthesis
p1 <-merge_gift_TM6 %>%
  ggplot(aes(x=type,y=Biosynthesis,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Type")

#Degradation
p2 <-merge_gift_TM6 %>%
  ggplot(aes(x=type,y=Degradation,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "type")

#Structure
p3 <-merge_gift_TM6 %>%
  ggplot(aes(x=type,y=Structure,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "type")

#Transport
p4 <-merge_gift_TM6 %>%
  ggplot(aes(x=type,y=Transport,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")

7.9.4.2 CI vs WC

#Merge the functional domains with the metadata
merge_gift_TM8 <- GIFTs_domains_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "Tube_code") %>% 
  inner_join(., sample_metadata_TM8, by="Tube_code")
#Biosynthesis
p1 <-merge_gift_TM8 %>%
  ggplot(aes(x=type,y=Biosynthesis,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Type")

#Degradation
p2 <-merge_gift_TM8 %>%
  ggplot(aes(x=type,y=Degradation,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 0.5, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "type")

#Structure
p3 <-merge_gift_TM8 %>%
  ggplot(aes(x=type,y=Structure,color=type,fill=type))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "type")

#Transport
p4 <-merge_gift_TM8 %>%
  ggplot(aes(x=species,y=Transport,color=species,fill=species))+
  geom_jitter(width = 0.2, size = 1.5, show.legend = FALSE)+ 
  geom_boxplot(alpha=0.2,outlier.shape = NA, width = 0.5, show.legend = FALSE, coef=0)+
  stat_compare_means() +
  theme(axis.text.x = element_text(vjust = 3, size=10),
        axis.text.y = element_text(size=10),
        axis.title=element_text(size=12,face="bold"),
        axis.text = element_text(face="bold", size=18),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        legend.text = element_text(size=10),
        legend.title = element_text(size=12),
        legend.position="none",
        legend.key.size = unit(1, 'cm'),
        strip.text.x = element_text(size = 12, color = "black", face = "bold"))+
  labs( x = "Species")